function obj = initPM100D(wavelen, varargin)

    p = inputParser;
    p.addRequired('wavelen',@(x)validateattributes(x,{'numeric'},{'>=',400,'<=',1100,'scalar','real'}));  % in unit of nm
    p.addParameter('SN','',@(x)ismember(x,{'P0014383'; 'P0021891'; 'P0016305';''}));
    % SN = 'P0014383'; % ZXD's device
    % SN = 'P0021891'; % SBW's device
    % SN = 'P0016305'; % YHY's device
    p.addParameter('avgN',10,@(x)validateattributes(x,{'numeric'},{'positive','scalar'}) );
    parse(p, wavelen, varargin{:});
    p = p.Results;
    obj = struct();
    
    % connect to device
    if strcmp(p.SN,'')
        obj.dev = instrfind('Type', 'visa-usb', 'ManufacturerID', '0x1313', 'ModelCode', '0x8078');   % Find all Thorlabs PM100D inst.
    else
        obj.dev = instrfind('Type', 'visa-usb', 'ManufacturerID', '0x1313', 'ModelCode', '0x8078' , 'SerialNumber', p.SN);
    end

    if isempty(obj.dev)
        if strcmp(p.SN,'')
            error('Cannot find PM100D device!! Please provide the SN printed on the device back panel (something like "P0014383").');
        else
            obj.dev = visa('NI', ['USB0::0x1313::0x8078::',p.SN,'::0::INSTR']);
        end
    else
        fclose(obj.dev);
        obj.dev = obj.dev(1);
    end
    fopen(obj.dev);
    fprintf(obj.dev, '*RST'); % reset
    pause(0.5);
    
    % record hardware info
    obj.IDN = query(obj.dev, '*IDN?','%s\n','%c\n');
    obj.sensor = query(obj.dev, 'system:sensor:idn?','%s\n','%c\n');

    % set wavelength
    fprintf(obj.dev, ['sense:correction:wavelength ',num2str(p.wavelen)]); 
    obj.wavelength = query(obj.dev, 'sense:correction:wavelength?','%s\n','%c\n');

    % set average number
    fprintf(obj.dev, ['sense:average ',num2str(p.avgN,'%u')]);
    obj.avgN = str2num(query(obj.dev, 'sense:average?','%s\n','%c\n')); %#ok<ST2NM>

    % configure measurement mode
    fprintf(obj.dev, 'configure:power');
    fprintf(obj.dev, 'measure:power');
    obj.measurementConfigure = query(obj.dev, 'configure?','%s\n','%c\n');

    % enable auto power range
    fprintf(obj.dev, 'sense:power:range:auto ON');
    obj.unit = query(obj.dev, 'sense:power:unit?','%s\n','%c\n');

    
    obj.readData = @()read(obj);
%     % start a new measurement & read data (recommended)
%     % PM100D starts to accumulate data when this command is sent, and return a result when the required number of samples is acquired.
%     % It takes 40~80ms when avgN = 100; 360ms for avgN = 1000; 3.4s for avgN = 10000; 3ms for avgN = 1; 7ms for avgN = 10.
%     tic;
%     query(obj.dev, 'read?','%s\n','%c\n')
%     toc;

%     % start a measurement, fetch data and then abort the measurement
%     % it takes about 3ms to fetch one result
%     fprintf(obj.dev, 'initiate');
%     for i = 1:5
%         tic;
%         query(obj.dev, 'fetch?','%s\n','%c\n')
%         toc;
%     end
%     fprintf(obj.dev, 'abort');
end

function P = read(obj)
    P = str2double(query(obj.dev, 'read?','%s\n','%c\n'));
end
